In [ ]:

    

import datetime


    

In [ ]:

    

# Load code to setup authorization for an IPython Notebook. Note that this is a temporary step
# that is required until the Earth Engine Python API is updated to include this logic.
%run 'authorize_earth_engine_in_notebook.ipynb'

# Initialize Earth Engine
# Note that we are calling a function defined in the previously run IPython Notebook, rather than
# the typical call to ee.Initialize()
ee_initialize()


    

In [ ]:

    

%run 'define_google_maps_interactive_widget.ipynb'


    

In [ ]:

    

# Define a comparison time interval.
testDateStart = datetime.datetime(2014, 1, 1)
testDateEnd = datetime.datetime(2015, 1, 1)

# Define a reference time interval.
referenceDateStart = datetime.datetime(2000, 1, 1)
referenceDateEnd = datetime.datetime(2014, 1, 1)

# Define date range within year that is used for comparison.
dayOfYearTarget = 10
dayOfYearInterval = 8 # interval length in days 

dayOfYearStart = dayOfYearTarget - 0.5 * dayOfYearInterval
dayOfYearEnd = dayOfYearTarget + 0.5 * dayOfYearInterval

base_collection_ndvi = ee.ImageCollection('MCD43A4_NDVI')
base_collection_ndsi = ee.ImageCollection('MCD43A4_NDSI')

# define the visualization parameters
vizNdvi = {
    'min': 0, 'max': 1,
    'palette': ','.join(
            ["FFFFFF","CE7E45","DF923D","F1B555","FCD163","99B718","74A901","66A000","529400",
             "3E8601","207401","056201","004C00","023B01","012E01","011D01","011301"])
  }
vizNdsi = {
    'min': -1, 'max': 1,
    'palette': 'FF0000,000000,0000FF'
  }
vizAnomaly = {
    'min':-0.4, 'max':0.4, 
    'palette': ','.join(
        ["87000A","7C3E28","EC712C","FABF45","FFFFFF","51FF78","3DCF4C","215229","0B260B"])
  }
thresholdAnomaly = 0.01 # suppresses the display of low values
vizLandcover = {
    'min':0, 'max':17,
    'palette': ','.join([
      'aec3d4', # water
      '152106', '225129', '369b47', '30eb5b', '387242', # forest
      '6a2325', 'c3aa69', 'b76031', 'd9903d', '91af40',  # shrub, grass, savanah
      '111149', # wetlands
      'cdb33b', # croplands
      'cc0013', # urban
      '33280d', # crop mosaic
      'd7cdcc', # snow and ice
      'f7e084', # barren
      '6f6f6f'  # tundra
    ])}

# create a land mask
landcover = (ee.Image('MCD12Q1/MCD12Q1_005_2009_01_01')
              .select('Land_Cover_Type_1'))
landMask = landcover.neq(0)

# create a agricultural area mask
agMask = landcover.eq(12)

# create a filter object that selects data based on the day-of-year
dayOfYearFilter = ee.Filter(
    ee.Filter.calendarRange(dayOfYearStart, dayOfYearEnd, 'day_of_year')
  )

# target collection
test_collection = (base_collection_ndvi
    .filterDate(testDateStart, testDateEnd)
    .filter(dayOfYearFilter))

# target collection - NDSI
test_collection_ndsi = (base_collection_ndsi
    .filterDate(testDateStart, testDateEnd)
    .filter(dayOfYearFilter))

# reference collection
reference_collection = (base_collection_ndvi
    .filterDate(referenceDateStart, referenceDateEnd)
    .filter(dayOfYearFilter))

anomaly = test_collection.median().subtract(reference_collection.median())

# Mask out non-agricultural areas, and snowy areas (NDSI>0)
#mask = agMask.and(test_collection_ndsi.median().lte(0))
mask = test_collection_ndsi.median().lte(0)

map = GoogleMapsWidget(lat=0, lng=0, zoom=2)
display(map)

map.addLayer(
    image=landcover,
    vis_params=vizLandcover,           
    name="landcover",
    visible=False)
map.addLayer(
    image=ee.Image(0),
    vis_params={'palette':"000000", 'opacity':0.5},           
    name="background",
    visible=True)
map.addLayer(
    image=reference_collection.median().mask(mask),
    vis_params=vizNdvi,           
    name="Median of Reference Collection",
    visible=False)
map.addLayer(
    test_collection.median().mask(mask),
    vis_params=vizNdvi,           
    name="Median of Test Collection",
    visible=False)
map.addLayer(
    test_collection.median().mask(mask),
    vis_params=vizNdsi,           
    name="Median of Test Collection NDSI",
    visible=False)
map.addLayer(
    anomaly.mask(anomaly.abs().gte(thresholdAnomaly).multiply(mask)),
    vis_params=vizAnomaly,           
    name="NDVI anomaly (test - reference)",
    visible=True)


    

In [ ]: